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AN  EXAMINATION  OF  STABILITY  CRITERIA  FOR  ITERATIVE  NUMERICAL  SCHEMES 
USED  IN  SOLVING  DIFFERENTIAL  EQUATIONS 


Katharine  Moore 


\  SUMMARY 

\  - 

jiiThe  stability  of  numerical  schemes  for  solving  algebraic  finite-difference 
equations  resulting  from  finite-difference  approximations  to  differential  equa¬ 
tions  is  discussed.  It  is  suggested  that  the  von  Neumann  method  together  with 
its  stability  criterion  provides  a  reasonably  simple  way  of  determining  stability. 
However,  there  are  limitations  in  its  applicability,  some  of  which  are  indicated. 
The  method  is  tested  in  two  examples  and  an  indication  is  given  of  how  best  to 
treat  first-  and  mixed-derivative  terms  occurring  in  differential  equations. 
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1  INTRODUCTION 

The  mathematical  modelling  of  practical  problems  often  involves  the  use  of 
differential  equations.  Very  few  of  these  equations  can  be  solved  analytically,  and  hence 
it  is  of  great  importance  to  develop  satisfactory  schemes  for  solving  them  numerically. 

One  requirement  of  any  satisfactory  numerical  scheme  is  that  it  should  be  stable.  How¬ 
ever  there  are  several  definitions  of  stability  in  the  literature,  all  leading  to  differ¬ 
ent  stability  criteria.  The  purpose  of  this  Report  is  to  discuss  some  of  these  defini¬ 
tions  and  criteria,  with  the  aid  of  examples.  The  examples  will  also  suggest  how  best  to 
treat  first-and  mixed-derivative  terms  when  they  arise  in  differential  equations. 

There  are  two  basic  steps  in  the  usual  finite-difference  methods  of  obtaining 
approximate  numerical  solutions  to  differential  equations.  First  a  finite-difference 
approximation  to  the  differential  equation  must  be  chosen,  the  result  of  this  being  a 
set  of  equations,  termed  the  finite-difference  equations.  The  second  step  is  to  solve 
these  equations,  thus  obtaining  an  approximation  to  the  solution  of  the  original 
differential  equation. 

Associated  with  these  two  steps  are  four  concepts:  consistency,  convergence, 
convergence  of  an  iterative  scheme,  and  stability.  These  will  be  explained: 

Consistency  If  the  finite-difference  formulation  is  equal  to  the  original  differential 
equation  plus  terras  which  tend  to  zero  as  the  grid  size  tends  to  zero,  then  the  finite- 
difference  approximation  is  consistent. 

Convergence  When  the  finite-difference  approximation  is  convergent,  the  difference 
between  the  discrete  approximation  (the  solution  of  the  finite-difference  equations)  and 
the  true  solution  (the  solution  of  the  differential  equation)  can  be  made  as  small  as 
desired,  by  choosing  a  sufficiently  small  grid. 

Convergence  of  an  iterative  scheme  If  the  finite-difference  equations  are  solved  by  an 
iterative  scheme,  this  scheme  is  convergent  if  and  only  if  the  sequence  of  approximations 
(the  (n  +  l)th  of  which  is  obtained  from  the  nth  by  the  iterative  scheme-)  converges  to 
the  solution  of  the  finite-difference  equations.  The  zeroth  approximation,  or  initial 
guess,  must  be  supplied. 

Stability  There  seems  to  be  much  disagreement  over  the  definition  of  this  term.  Some 

authors  merely  require  that  all  errors  should  eventually  be  damped  out'.  Others  appear 

.  2  3 

to  relate  stability  to  the  growth  of  rounding  errors  *  . 

For  initial-value  problems  Richtmyer  and  Morton*  give  several  definitions  of 
stability.  These  do  not  mention  errors,  as  they  are  given  in  terms  of  the  operators  by 
which  the  solution  at  time  t  +  At  may  be  obtained  from  the  solution  at  time  t  ,  where 
'time'  is  the  coordinate  in  the  marching  direction,  and  At  the  time  step  taken.  These 
'operator'  definitions  can  be  extende..  to  iterative  schemes  for  solving  sets  of  equations 
arising  from  finite-difference  approximations  to  elliptic  partial  differential  equations, 
if  the  time-dependent  analogy  of  Jameson  (see  for  instance,  Ref  4)  is  used.  In  this 
analogy  an  artificial  time  coordinate,  t  ,  with  step  length  At  ,  in  introduced,  and  <j>n  . 
the  approximation  to  the  solution  after  the  nth  step,  is  regarded  as  the  solution  at  time 
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The  last  two  ideas  depend  only  on  the  finite-difference  equations  and  the  manner  in 
whicn  they  are  solved,  no  reference  being  made  to  the  original  differential  equation. 

They  ensure  that  the  solution  of  the  finite-difference  equations  can  be  obtained. 

There  are  many  theorems  which  show  that  convergence  (and  convergence  of  the 
iterative  scheme  if  appropriate)  will  be  obtained  for  various  classes  of  differential 
equations  with  fairly  general  initial  and  boundary  conditions,  if  the  finite-difference 
equations  are  consistent  and  the  scheme  adopted  for  their  solution  is  stable  in  a  suitable 
sense  (see,  for  instance,  Ref  1).  However,  for  nonlinear  second-order  partial  differen¬ 
tial  equations  there  are  no  such  results  except  in  a  few  special  cases.  Nevertheless 
for  such  equations  it  is  widely  accepted  (and  will  be  assumed  here  from  now  on)  that 
stability  (in  some  suitable  sense)  and  consistency  do  imply  convergence  (and  convergence 
of  the  iterative  scheme  if  appropriate) . 

It  is  usually  straightforward  to  show  consistency  -  in  practice  this  is  done  in  the 
formulation  of  the  finite-difference  approximation.  In  section  2  the  discussion  on 
stability  will  be  continued  and  it  is  suggested  that  the  von  Neumann  criterion  is  a 
reasonably  simple  one  which  is  adequate  in  many  cases.  Attention  is  also  drawn  to  some 
of  its  limitations.  In  section  3  some  of  the  ideas  introduced  in  section  2  will  be 
illustrated  with  specific  examples,  and  it  is  shown  how  firsthand  mixed-derivative 
terms  might  best  b>.  treated. 

2  STABILITY  OF  ITERATIVE  SCHEMES 

For  many  types  of  differential  equations  (for  example,  wnen  the  coefficients  of 
the  highest  derivative  terms  are  functions  of  the  solution)  there  are  no  rigorous  theories 
concerning  the  stability  of  numerical  schemes  for  their  solution.  The  usual  approach  is 
to  do  a  local  stability  analysis,  and  to  hope  that  if  the  scheme  is  everywhere  locally 
stable  it  will  be  globally  stable.  Although  this  is  not  always  true,  practical  experience 
suggests  the  correspondence  is  close,  probably  because  instabilities  usually  arise 
locally' . 

There  are  two  methods  commonly  used  for  examining  the  notion  of  stability  of  a 
.  3 

finite-difference  scheme  .  The  first  one,  termed  the  von  Neumann  method,  will  be 
examined  in  section  2.1.  The  second  one,  termed  the  matrix  method,  will  be  discussed  in 
section  2.4 

2. 1  The  von  Neumann  stability  method 

The  differential  equation  is  taken  to  have  constant  coefficients,  and  the  problem 
is  assumed  to  be  an  initial  value  one,  the  only  permitted  boundary  conditions  being  ones 
which  can  be  replaced  by  periodicity  conditions.  The  'amplification  factors',  X  , 
can  then  be  determined  as  follows: 

(i)  Substitute,  into  the  finite-difference  equations  used  to  solve  the  differential 

equation, 

un(x)  -  Xn^  exp  ( ik  ,  x) 


where  k  is  a  real  vector  chosen  to  satisfy  the  periodicity  conditions,  and 
un(x)  is  the  vector  of  unknowns  at  position  x  after  the  nth  step. 

(ii)  Determine  the  amplification  factors  from  the  condition  for  the  existence  of 
a  non-zero  solution  for  Uq  . 

von  Neumann's  criterion  states  that  a  necessary  condition  for  stability  (as  defined  in 
Richtmyer  arid  Morton*,  Chapter  3)  is  that  all  the  amplification  factors,  X  ,  must 
satisfy 

|  X |  <  1  +  0(h)  for  0  <  h  <  t 

where  t  is  some  upper  bound  on  h  ,  and  h  is  the  step  length  in  the  marching 
direction.  This  stability  criterion  can  be  extended  to  the  Iterative  methods  for 
solving  elliptic  partial  differential  equations  with  periodic  boundary  conditions,  by 
use  of  Jameson's  time-dependent  analogy  (see,  for  instance,  Ref  4).  It  then  becomes 
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The  application  of  the  von  Neumann  criterion  will  be  illustrated  by  the  following 
example*,  in  which  the  boundary  conditions  are  periodic: 

2 

solve  "  a  “•?  on  0  <  x  <  1  and  0  <  t 

9t  3x2 

with  a  a  constant,  u(0,t)  ■  u(l,t)  -  0  and  u(x,0)  ■  F(x)  .  Let  u"  be  the  finite- 
difference  solution  at  x  *»  jAx  and  t  »  nAt  ,  take  a  central-difference  approximation 

2  2  1^4.  j 

for  3  u/Ex  ,  and  use  an  explicit  scheme  (that  is  one  in  whi  .h  one  unknown  u.  is 

n  ^ 

expressed  in  terms  of  the  known  u.s).  The  finite-difference  equations  may  then  be 

written 


n+1  n 
u.  -  u. 

-J _ 1 

At 


-A 

3 


(Ax) 


I  n  0  n  .  n  \ 
u.  ,  -  2u .  +  u.  ,) 
2  \  J+l  J  1-1/ 


giving  the  equations,  for  the  amplification  factors,  X  , 

ikAx 


(X  ~  '.) 
At 


o  ikAx  -ikAx 
- r  e  +  e  -  2 

(Ax)2  L  J 


!  I 

i  1 

:i 


where  k  *  2wn  ,  with  n  any  integer,  to  satisfy  the  periodicity  condition.  Hence 


ao 

© 


1 - „  2(1  -  cos  kAx] 

(Ax) 2 


and  the  von  Neumann  criterion  gives 


2oAt 

(Ax)2 


<  1 


as  a  necessary  condition  for  stability. 
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Under  some  circumstances  the  von  Neumann  condition  is  sufficient  as  well  as  neces- 
s ary  for  stability.  In  particular,  for  two-level  schemes  (that  is  schemes  in  which  the 
finite-difference  equations  relate  values  of  the  dependent  variables  at  the  (n  +  l)th 
step  to  values  at  the  nth  step,  values  of  the  dependent  variables  at  earlier  steps  not 
occurring  in  the  equations)  with  one  dependent  variable  and  any  number  of  independent 
variables,  the  von  Neumann  criterion  is  sufficient  as  well  as  necessary  for  stability. 

For  U.e  example  given  above  equation  (2-1)  is  thus  a  necessary  and  sufficient  condition 
to  ensure  stability, 

In  cases  where  the  von  Neumann  criterion  is  necessary  but  not  sufficient  for 
stability,  further  criteria  which  ensure  stability  can  often  be  found1 . 

Although  it  might  appear  here  that  the  von  Neumann  method  has  limited  applicability, 
particularly  because  of  the  restriction  to  periodic  boundary  conditions,  in  fact  it  has 
much  wider  practical  application.  In  section  2.4.2  the  Godunov-Ryabenkii  criterion  is 
described.  This  in  effect  provides  an  extension  to  the  von  Neumann  treatment  to  cover 
consideration  of  arbitrary  boundary  conditions,  although  in  practice  the  extension  is 
often  difficult  to  apply. 

2.2  Matrix  formulation  of  the  finite-difference  problem 

The  problem  of  finding  the  solution,  $  ,  of  a  set  of  finite-difference  equations 
approximating  a  differential  equation  (ordinary  or  partial)  can  be  expressed  in  the 
form: 

find  <j>  satisfying  A<J>  «■  B  (2-2) 

where  B  is  independent  of  $  ,  but  A  may  depend  on  4  .  As  an  example  consider 
the  numerical  solution  of 

a4  +  2b4  +  C4  ■  f 

xx  xy  yy 


on  a  unit  square,  with  Ax  ■  Ay  “  1/N  ,  a,  b,  c,  f  constants  and  values  of  4  given  on 
the  boundary.  Let  the  subscripts,  i  and  j  ,  refer  to  the  coordinate  directions,  x 
and  y  ,  respectively.  Take  the  usual  central-difference  representations  of  4>xx  and 
♦  and  assume  that  if  *  were  known,  4^  would  be  approximated  by 


*  (iAx.jAy) 


4AxAy  ^*i+l,j+l  ^i+l,j-l 


Vi.j+i  1-1,  j- 


(2-3) 


where  4^  ■  4(iAx,jAy)  . 

The  equation  (2-2)  takes  the  form 


a 

3 
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where  each  Yj  is  a  vector  of  length  (N  -  1).  The  ith  component  of  the  vector  is 

$  [i  +  (j  -  1)(N  -  1)J  and  is  the  estimate  which  is  obtained  for  from  the  finite- 

difference  approximation.  B  is  a  vector  of  length  (N  -  1 ) ~  with  elements  equal  to 
2 

f/N  plus  (possibly)  some  contribution  from  the  boundary.  E  and  F  are  matrices  of 
order  (N  -  1)  x  (N  -  1)  with 


and 


(2-5) 


T  T 

F  denotes  the  transpose  of  F  ,  that  is  F^  ■  Fji  ’  The  Partlcu^ar  form  of  the 

matrix,  A  , ^ which  is  of  order  (N  -  1)^  *  (N  -  l)^j  is  determined  by  the  finite- 
difference  approximations  used. 

The  matrix,  A  ,  has  several  features  of  interest.  It  is  sparse,  which  means  that 
most  of  its  entries  are  zero  -  there  is  no  precise  definition  of  sparse,  but  a  good 
guide  seems  to  be  that  a  matrix  is  sparse  if  more  than  90Z  of  its  entries  are  zero.  It 
is  block  tri-diagonal,  and  each  non-zero  block  is  tri-diagonal.  Features  of  this  nature 
are  generally  observed  with  finite-difference  equations  obtained  from  partial  differential 
equations. 

2. 3  Methods  of  solution 

Consider  equation  (2-2)  again.  For  any  given  A  this  set  of  equations  can  be 
solved  directly  by  Gaussian  elimination.  Frequently  some  form  of  reordering  of  A  before 
carrying  out  the  elimination  will  give  a  much  faster  scheme.  This  sort  of  scheme  will 
be  most  suitable  when  A  is  independent  of  $  (as  will  occur  if  the  differential 
equation  is  linear),  or  if  A  is  triangular  (either  upper  or  lower)  with  all  elements 
on  the  leading  diagonal  independent  of  $  (as  might  occur  if  the  differential  equation 
is  hyperbolic  or  parabolic). 

However,  the  equations  can  also  be  solved  iteratively  by  writing  A  in  the  form 


■J  t JaLiiSlM  l-  Jt. . . :  i  .A.,  !,  ‘ir  T^i^iljiV  T  -  :  ..iu^Uakak^ 


8 


and  solving 


Gw 


n+1 


B  +  Cwn 


(2-6) 


where  w^  is  given,  G  is  some  easily  invertible  matrix  and  any  appearing  in  G 
and  C  is  replaced  by  the  estimate  w”  of  ^  .  If  the  scheme  is  stable  in  the 
sense  of  all  errors  decaying  suitably  then  wn  will  tend  to  $  as  n  ■+  «  .  An  itera¬ 
tive  scheme  of  this  sort  will  obviously  be  most  suitable  when  A  is  dependent  on  $  , 
and  equation  (2-2)  cannot  be  solved  directly,  as  might  occur  if  the  differential  equation 
is  elliptic  and  'quasi-linear' .  A  quasi-linear  differential  equation  is  one  in  which  the 
highest  derivatives  appear  linearly,  but  the  coefficients  of  the  highest  derivatives  are 
functions  of  the  dependent  variable  and  lower-order  derivatives  of  it. 

It  is  of  interest  to  note  that  approximating  some  second-order  parabolic  and 
first-order  hyperbolic  quasi-linear  partial  differential  equations  by  certain  'marching' 
finite-difference  schemes  can  also  give  rise  to  equations  of  the  form  (2-'),  where  wn 
is  now  the  approximation  to  the  solution  at  the  nth  step  in  the  marching  direction. 

This  is  illustrated  in  the  example  described  in  subsection  2.4.2.  Equation  (2-2)  now 
takes  the  form 


2.4  Matrix  method  of  determining  stability 

2.4.1  Preliminaries 

It  was  indicated  at  the  beginning  of  this  section,  that  there  are  very  few 
theories  concerning  stability,  when  the  matrix,  A  ,  of  equation  (2-2)  depends  on  the 
solution  of  the  finite-difference  equations.  A  local  stability  analysis  is  usually 
carried  out  in  such  cases.  Hence,  from  now  on,  the  matrices  A  of  equation  (2-2),  and 
G  and  C  of  equation  (2-6)  will  be  assumed  to  be  independent  of  the  solution,  ,  of 
the  finite  difference  equations. 

If  the  'error',  e11  ,  is  now  defined  by 


it  can  easily  be  seen  that 


Ge 


n+1 


Ge 


(2-7) 


o 

00 


i 


i 


and  that  w  will  tend  to  $  if  and  only  if  e11  tends  to  zero,  as  n  tends  to 
infinity.  Thus,  it  might  be  expected  that,  in  principle  at  least,  stability  could  be 
discussed  in  terms  of  the  properties  of  the  matrix  G^C  . 


jJaimifeLBjyillf.  -l>..f.<l,  ,i^i.  *  ...tu..  i:.  «  ..i: _ i  *Ml  .  ,.r  .  -  -.l..  J'li.ji.  K  rH.-.r*iiUJ;i  1.  ml.  *i.  ..i.  ..  j  .1  ...A . ... ..t*  L  iA^,,jlliid|4litife.i..i.4z  ,1 
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Before  diecueeing  some  possible  stability  criteria  it  is  necessary  to  define  a  few 
terms  and  some  notation: 


xi*f  »  where  x  is  a  vector  of  length  n  ,  and  *  denotes  complex  conjugate. 


Hermitian  A  square  matrix,  A  ,  is  hermitian  if  A^  ■  A*j  . 

Diagonally  dominant  An  n  *  n  complex  matrix,  A  ,  iB  diagonally  dominant  if 


i Aii I  > 


t  'V 


for  all  i  . 


Negative  definite  An  hermitian  n  *  n  matrix,  A  ,  is  negative  definite  if,  for  all 
vectors  x  not  identically  zero,  and  of  length  n 


it** 


<  0  . 


i-1  l-l 


Eigenvalue  and  corresponding  eigenvector  X  is  said  to  be  an  eigenvalue  (also  termed 
latent  root  or  characteristic  root)  of  an  n  *  n  matrix  A  ,  if  there  exists  a  non-zero 
vector  e  of  length  n  such  that 

Ae  -  Xe  , 


The  vector  e  ,  which  is  determined  to  at  most  an  arbitrary  multiplicative  constant,  is 
termed  the  eigenvector  corresponding  to  the  eigenvalue,  X  . 

Orthogonal  Two  vectors,  x  and  ^  ,  each  of  length  n  are  orthogonal  if 


x  •  l 


tv, 


Spectral  radius,  p(A)  The  spectral  radius  of  an  n  *  n  complex  matrix  is  p(A)  where 

P (A)  -  max  ( X . | 

!  <  i  <  n 

and  {X.}  is  the  set  of  eigenvalues  of  A  . 

Spectral  norm,  |a|  The  spectral  norm  of  an  n  *  n  complex  matrix  is  II aH 


where  U  All  -  max 

all  x^O 


and  x  is  a  vector  of  length  n 
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2.4.2  Stability  criteria 

Three  possible  criteria  will  1  .  described  in  this  section  and  a  discussion  of  them 
follows  in  section  2.4.3. 

Clearly,  if,  in  equation  (2-7)  |en|  +0  m  n  +  » 

p(g"*C)  <  1  .  (2-8) 

In  equation  (2-6)  this  will  ensure  that  errors  eventually  decay,  so  stability  is  achieved 
in  one  of  the  senses  of  section  1 .  For  the  rest  of  this  Report  condition  (2-8)  will  be 
termed  the  first  stability  criterion.  This  criterion  seems  t.o  be  widely  used  -  see, 
for  instance,  Refs  5  and  6. 

The  most  obvious  shortcoming  of  this  criterion  is  that  for  parabolic  and  hyperbolic 
equations  it  will  exclude,  the  possibility  of  solutions  which  are  growing  exponentially 
in  the  marching  direction.  To  include  these,  the  stability  criterion  should  be  relaxed 
to 

p(G~*C)  <  1  +  0(h) 


where  h  is  the  step  length  in  the  marching  direction.  However,  stability,  in  any  of 
the  senses  given  in  section  1,  is  no  longer  necessarily  achieved. 

Examples* show  further  unsatisfactory  features  of  the  stability  criterion 
(2-8).  This  will  be  illustrated  here  with  the  example': 

solve  -5^  +  a  -5^  "0  on  0  <  x  <  1  and  0  <  t  , 


With  a  a  positive  constant,  u(0,t)  -  0  and  u(x,0)  ■  F(x)  .  The  true  solution 
is  u  «  F(x  -  at)  when  x  >  at  ,  and  zero  everywhere  else.  Take  Ax  *  1 /N  ,  represent 
3u/3x  by  a  backward  difference,  and  use  an  explicit  scheme.  This  gives 


n+1  n 
u.  -  u 


J _ 1  .  JL  (un  .  un  \ 

At  Ax  \  j  j-I/ 


where  u?  is  the  calculated  value  of  u  at  x  ■  jAx  and  t  *  nAt  .  Let  r  »  Ata/Ax  , 
then 

uf1  -  (1  -  r)u?  +  ru?  ,  (2-9) 

J  J  J 


giving  the  matrix  G  as  the  identity  and  C  as  the  N  *  N  matrix 


oo 

o 


For  fixed  i  it  can  be  shown  that  |e?|  has  a  maximum  at  q  ■  2i  -  2  ,  and  that  this 
maximum  is  an  increasing  function  of  i  .  Hence 


max 

i.q 


m 


’l\N"'/3\N~'  (2N  -  2)1  . 

JJ  \2  )  (N  -  1 )  !  (N  -  ITT  1 


which  is  asymptotic  to,  for  large  N  , 
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by  Stirling's  formula.  This  shows  that  stability  is  not  achieved  in  the  sense  of  the 
effect  of  a  rounding  error  being  bounded  as  N  increases,  nor  is  it  achieved  in  the 
sense  of  any  of  the  definitions  given  in  Richtrayer  and  Morton' .  In  numerical  work  the 
scheme  will  become  less  and  less  satisfactory  as  Ax  ■  1 /N  becomes  smaller  and  smaller, 
while  r  is  kept  fixed.  However,  for  small  enough  N  ,  it  may  appear  satisfactory 
because  max  |ei^|/e  will  not  be  very  large,  and  errors  increasing  only  slightly  before 

i.q 

decaying  may  be  acceptable. 

Mathematically,  the  difficulties  of  errors  becoming  arbitrarily  large  before  they 
decay  away,  as  the  grid  size  decreases,  seem  to  arise  either  through  there  being  fewer 
essentially  distinct  eigenvectors  than  eigenvalues  of  G  'c  (as  in  the  example  just 
discussed),  or  through  the  eigenvectors  being  nowhere  near  orthogonal^ ' 7 . 

To  prevent  components  of  any  introduced  error  becoming  arbitrarily  large  as  the 
grid  size  is  decreased,  one  might,  as  suggested  in  Richtmyer  and  Morton' ,  impose 


TITO?-*  ?fPWSB?  » If •FT ’’Tpsr  tTwr, •rir  :5 
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some  constant, 


(2-10) 


independent  of  the  integer,  v  ,  and. the  grid  size  for  all  sufficiently  small  grid  sizes, 
v  is  any  integer  satisfying  0<hv  <H  ,  where  H  is  some  given  constant,  and  h  is 
the  step  length  in  the  marching  direction,  which  satisfies  0  <  h  <  t  ,  x  being  some 
suitable  upper  bound.  This  condition  can  be  extended  to  iterative  methods  for  solving 
elliptic  partial  differential  equations,  by  use  of  Jameson' a  time-dependent  analogy  (see, 
for  instance.  Ref  4).  It  is  the  same  as  that  given  in  equation  (2-10),  except  that 
v  is  now  any  positive  integer. 

Satisfying  this  constraint,  which  will  be  termed  the  second  stability  criterion 
for  the  rest  of  this  Report,  seems  to  ensure  stability  in  all  the  senses  given  in 
section  1 .  However,  for  this  constraint  to  be  satisfactory  in  practice,  the  constant, 
k  ,  must  not  be  too  large,  and  the  grid  size  must  not  be  too  small. 

In  the  example  discussed  earlier  in  this  section  this  stability  criterion  gives 
0  <  r  <  1  ,  which  is  the  well-known  Courant-Friedrichs-Lewy  criterion  for  the  stability 
of  the  iterative  scheme. 


If  there  is  only  one  space-like  independent  variable  it  can  be  shown  that,  if  the 
second  stability  criterion  is  to  be  satisfied,  it  is  necessary  that  the  Godunov-Ryabenkii 
criterion  be  satisfied.  A  statement  of  the  Godunov-Ryabenkii  criterion  and  a  proof  of 
its  necessity  can  be  found  in  Richtmyer  and  Morton*.  It  is  in  two  parts,  one  part 
being  a  condition  that  looks  very  like  a  von  Neumann  criterion.  Amplification  factors, 

X  ,  (local  ones  if  necessary)  are  found  by  exactly  the  same  procedure  as  described  in 


section  2.1.  These  must  satisfy 

lim  | A  |  <  1 

h-H) 

where  h  is  the  step  length  in  the  marching  direction. 

This  bound  on  the  sizes  of  the  amplification  factors  is  very  similar  to,  but  weaker 
than,  the  bound  imposed  on  th^m  in  section  2,i.  By  use  of  the  time-dependent  analogy 
(see,  for  instance,  Ref  4)  this  ’.icomes 

W  <  1 

for  iterative  methods  used  in  solving  elliptic  partial  differential  equations.  This 
bound  on  the  sizes  of  the  amplification  factors  is  the  same  as  that  imposed  on  them  in 
section  2.1. 

In  the  example  discussed  earlier,  equation  (2-9)  can  be  used  to  show  that  the 
amplification  factors,  X  ,  have  the  form 


0  <  r  <  1  . 


f  t  The  second  part  of  the  Godunov-Ryabenkii  criterion  is  concerned  with  the  boundary 

f  ¥  conditions,  each  boundary  condition  being  treated  in  turn.  The  implementation  of  this 

f  f  part  of  the  criterion  for  the  general  case  is  discussed  in  Richtmyer  and  Morton  ,  How- 

[  |  ever  the  discussion  can  be  much  simplified  if  the  following  restrictions  are  imposed: 

t  ? 

I  (i)  The  differential  equation  has  constant  coefficients. 


(ii)  The  differential  equation  has  a  scalar  dependent  variable. 

(iii)  The  highest  space  derivative  is  at  most  second  order. 

(iv)  A  two-level  difference  scheme  (this  term  is  defined  in  section  2.1)  is 
adopted. 

(v)  There  is  at  most  one  boundary  condition  at  each  boundary. 

Attention  here  will  thus  be  confined  to  this  class  of  problem. 

There  must  be  two  boundary  conditions  for  a  problem  in  which  the  highest  space 
derivative  is  of  order  two  and  one  for  a  problem  in  which  the  highest  space  derivative 
is  of  order  one  -  it  is  taken  that  this  is  true.  For  clarity  it  will  be  further  assumed 
that  the  problem  has  been  formulated  in  terms  of  the  equation  satisfied  by  the  errors 
and  the  boundary  conditions  on  them.  Errors  must  satisfy  homogeneous  boundary  conditions 
( ie  errors  do  not  contribute  anything  to  the  boundary  conditions) ,  because  the  Godunov- 
Ryabenkii  criterion  is  only  applicable  to  problems  with  linear  boundary  conditions. 

If  a  particular  boundary  condition  on  the  ^oendent  variable  is  homogeneous,  errors 
also  satisfy  this  boundary  condition. 

If,  when  carrying  out  a  von  Neumann  stability  analysis  on  the  equation  describing 
the  behaviour  of  the  errors,  exp(ikAx)  is  replaced  by  y  ,  an  equation  relating  X  and 
V  results.  (To  carry  out  the  von  Neumann  stability  analysis  this  equation  must  now  be 
solved  with  y  =  exp(ikAx)  ,  for  all  possible  values  of  kAx .)  If  the  substitutions 


en(j Ax)  =  Xnyje0  (2-11) 

where  en  is  the  error  at  the  nth  step  after  it  is  introduced  and  e^  is  a  non-zero 
constant,  are  now  made  into  the  finite-difference  equation  modelling  the  boundary 
condition  on  the  error,  at  the  boundary  under  consideration,  another  equation  for  X  and 
y  results.  These  two  equations  can  be  solved  simultaneously  for  X  and  y  ;  |y|  need 


not  necessarily  equal  unity. 

If  the  value  of  u  ,  the  dependent  variable,  is  specified  at  the  boundary,  the  error 
must  always  be  zero  at  the  boundary,  and  so  y  must  be  zero,  because  e^  cannot  be  zero. 
This  is  always  stable.  In  cases  when  y  is  non-zero,  e  ,  given  in  equation  (2-11)  is 
only  acceptable  (physically)  if  it  decays  away  from  the  boundary  under  consideration, 
into  the  interior  of  the  space*.  Hence  at  a  lower  boundary  only  solutions  for  X  and  u 


. . . gjJlAiAfab*. 


with  |y|  <  1  are  acceptable,  while  at  an  upper  boundary  only  solution i  with  |p|  1 

are  acceptable.  The  boundary  condition  under  consideration  is  said  to  be  stable  if 
y  ,  corresponding  to  an  acceptable  u  ,  satisfies 

lim  |X|  <  1 

h+0 

for  marching  problems,  with  h  the  step  length  in  the  marching  direction,  and 

N  <  i 

for  elliptic  problems.  The  second  part  of  the  Godunov-Ryabenkii  criterion  is  satisfied 
if  all  boundary  conditions  are  stable. 

It  is  often  rather  more  difficult  to  ascertain  whether  this  part  of  the  Godunov- 
Ryabenkii  criterion  is  satisfied,  than  it  is  to  ascertain  whether  the  'von  Neumann- like' 
part  of  the  criterion  is  satisfied.  However,  in  the  example  earlier  in  this  section, 
the  boundary  condition  at  x  •  0  cannot  lead  to  any  instabilities,  because  the  dependent 
variable  is  zero  there. 

As  has  already  been  commented  the  second  stability  criterion  might  not  be 
sufficiently  strong.  To  avoid  the  possibility  that  the  constant,  <  ,  of  equation  (2-10) 
might  be  too  large  in  practice,  one  might  impose  what  will,  for  the  rest  of  this  Reoort, 
be  termed  the  third  stability  criterion, 

II  G~ 1 C  H  <  1  +  0(h)  (2-12) 

where  h  is  the  step  length  in  the  marching  direct’’ on  (which  is  to  be  regarded  as  zero 
when  there  is  no  marching  direction).  This,  however,  seems  over-restrictive,  because 
it  very  rarely  matters  if  an  error  does  grow  slightly  before  decaying. 

2.4.3  Discussion  of  criteria 

It  should  be  noted  that,  if  the  third  stability  criterion  (2-12)  is  satisfied, 
then  the  other  two  are  necessarily  satisfied  because  (2-12)  (2-10)  (2-8).  On  the 

other  hand,  if  the  eigenvectors  of  G  'c  are  mutually  orthogonal  and  there  are  as  many 
eigenvectors  as  there  are  rows  of  G  C  ,  then,  if  the  first  stability  criterion  (2-8) 
is  satisfied,  the  other  two  are  necessarily  satisfied.  This  is  shown  by  proving 

(2-8)  *♦  (2-10)  "*•  (2-12) 

under  such  circumstances. 

For  practical  purposes  what  is  desired  is  a  necessary  and  sufficient  condition  for 
stability  (in  some  suitable  sense)  which  is  easily  applied  to  any  numerical  scheme. 

For  most  problems  it  is  very  difficult  to  find  p(G  *C)  or  II  (G  *0^11  for  all  permissible 
integers,  v  .  For  many  problems  with  one  space  dimension  it  is  also  difficult  to  apply 
the  part  of  the  Godunov-Ryabenkii  criterion  involving  the  boundary  conditions.  (The 
Godunov-Ryabenkii  criterion  is  only  applicable  to  problems  with  one  space-like  independent 
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variable.)  Only  an  analysis  of  the  form  described  in  section  2.1,  leading  to  a 
'von  Neumann-like'  stability  criterion  is  usually  fairly  straightforwtrd.  It  has  so 
far  only  been  indicated  that  'von  Neumann-like'  criteria  are  applicable  under  rather 
restrictive  conditions  (see  sections  2.1  and  2.4.2).  However,  it  is  widely  accepted* 
that  criteria  of  this  form  are  applicable  to  many  problems  not  included  in  the  categories 
described  in  sections  2.1  and  2.4.2.  Intuitively  this  'feels  right'  because,  away  from 
the  boundaries,  such  criteria  give  necessary  conditions  for  local  stability,  and  the 
correspondence  between  local  stability  and  global  stability  is  usually  close*.  However, 
as  the  Godunov-Ryabenkii  criterion  indicates,  even  if,  in  a  von  Neumann-like  analysis 
the  amplification  factors  all  have  moduli  less  than  unity,  the  scheme  may  still  be 
unsatisfactory  in  practice,  because  of  the  boundary  conditions.  This  will  be  illustrated 
in  section  3.2.  In  contrast,  in  some  cases  where  the  amplification  factors  do  not  obey 
a  von  Neumann-like  criterion  the  numerical  scheme  may  appear  to  be  satisfactory.  This 
will  arise  if  the  grid  is  sufficiently  coarse  and  p(G  *C)  <1  .  An  example  of  this 
will  be  given  in  section  3.1.  In  these  circumstances  the  scheme  will  become  less  and 
less  satisfactory  as  the  grid  is  refined. 

As  the  only  useful  criterion  that  has  been  suggested  (ie  a  von  Neumann-like 
criterion)  is  a  necessary  condition  for  stability,  it  might  be  of  value  to  see  if  useful 
criteria  which  are  sufficient  for  stability  can  be  found.  As  mentioned  in  section  2.1, 
conditions  which,  with  the  von  Neumann  criterion,  are  sufficient  for  stability  for 
schemes  approximating  pure  initial  value  problems  can  be  found  in  Richtmyer  and  Morton*. 
For  problems  involving  boundary  conditions  and  variable  coefficients  the  most  useful 
method  for  finding  sufficient  conditions  for  stability,  and  incidentally  for  indicating 
suitable  numerical  schemes,  is  probably  the  so-called  'energy  method'*.  However  the 
application  of  the  method  usually  seems  to  require  much  algebra,  and  often  leads  to 
very  complicated  conditions  which  are  far  from  necessary. 

A  suitable  practical  approach  to  obtaining  stability  might  thus  be  to  ensure  that 
the  von  Neumann-like  criterion  of  section  2.4.2  is  satisfied  locally  everywhere. 

3  EXAMPLES 

Here  the  ideas  of  section  2  will  be  illustrated  with  two  examples.  In  section  3.1 
the  representation  of  a  mixed  second  derivative  in  a  second-order  elliptic  differential 
equation  will  oe  discussed.  In  section  3.2  a  second-order  ordinary  differential  equation, 
in  which  occur  terms  involving  the  first  derivative,  will  be  considered.  Suitable 
iterative  schemes  for  solving  the  finite-difference  equations  will  be  suggested. 

3. 1  Example  1  -  mixed  second  derivative 

Consider  the  equation 

a4>  +  2b|t  +  ct  +  d<l  »  0 

xx  xy  yy  x 

2 

on  a  unit  square  with  $>  zero  everywhere  on  the  boundary  and  ac  >  b  (ie  the  equation 
is  elliptic).  Without  loss  of  generality  take  a  to  be  positive.  The  only  solution  of 
the  equation  is  $  identically  zero  everywhere.  Let  the  subscripts,  1  and  j  ,  refer 
to  the  coordinate  directions,  x  and  y  ,  respectively,  and  let  Ax  ■  Ay  »  1/N  .  Take 
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the  usual  central  difference  representations  of  *  ,  *  and  ♦  .  One  possible 

xx  yy  ^  # 

representation  of  the  crosa-derivativa  was  given  in  equation  (2-3) .  This  representation 
has  a  truncation  error  of 


4-  (Ax)2*  +  (Ay)2*  +  higher  order  terms. 

6  J_  xxxy  J  xyyyj 


Mitchell0  suggests 


*  (iAx,jAy)  —  ..  . 
xy  ,J  J/  2AxAy 


[*  -  *  -  * .  .  +  * . 

i+l,j+l  i,j+l  i+l. J  i.J 

+  *,  .  -  *.  ,  .  -  *.  .  ,  +  *.  .  .  .1 
l.J  1-1.  J  l.J“l  T--1»J"1J 


if  b  >  0  (3-la) 


with  *„  -  *(iAx,jAy),  which  has  a  truncation  error 


—  —  —  $  +  ■■  t  +  higher  order  terms 

6  xxxy  A  xxyy  6  xyyy 


and 


*  (iAx.jAy)  -  -  — 

xy  ,J  J  2AxAy 


A  -  -<t)  +  <J> 

Li.j+i  Vi,  j+i  vij  l-i  ,j 

+  *.  .  .-*..-  ♦.  .  ..+*..  ."1  if  b  <  0  (3-lb) 

i+l. J  1J  i+l, 3-1  i.J-lJ 


3  •  . 

with  a  similar  truncation  error.  Mitchell  suggests  this  scheme  so  that  the  coefficients 

of  all  the  unknowns  in  the  resulting  finite-difference  equations  will  be  positive.  He 

therefore  requires  |b|  <  min(a,c)  ,  in  the  case  of  zero  d  .  While  this  sort  of  condi- 

3  8 

tion  is  necessary  for  some  of  the  matrix  theories  of  stability  that  have  been  developed  ’ 
it  is  not  usually  a  necessary  condition  for  stability  and  is  of  no  use  in  the  situation 
b  >  min(a,c)  .  It  is  thus  worthwhile  exploring  alternative  representations  of  the 
cross-derivative . 

The  third  and  last  representation  of  the  cross-derivative  to  be  considered  here 
will  be  termed  the  'inverse  Mitchell'  scheme,  because  it  is  the  same  as  Mitchell's 
scheme  (given  above  in  equations  (3-1)),  ^except  that  the  condition  on  b  is  inverted. 

It  is 


VlAx.jAy)  3  2^a7 


i » j + 1  *i-l»j  +  l  *i.j  +  *i-l»j 


+  *i+l»j  "  *i»j  *i+l.j"l  +  *i.j“l] 


if  b  >  0  (3-2a) 


o 

00 


if  b >0  (3-2b) 


*xy^lAX,jAy  2AxAy  Li+J»i+1  ~  i*j  +  1  "  i+l»J  i»j 


i,J  i-l.J  i.J->  i-I.J-'J 


These  three  possible  representations  for  the  cross-derivative  all  lead  to  finite- 
difference  equations  of  the  form 

A$  «  0  (3-3) 

where,  as  in  section  2.2,  <J>  ti  +  (j  -  1)(N  -  1)]  is  the  estimate  of  $  at  x  ■  iAx  , 

y  -  j  Ay  which  is  obtained  from  solving  the  finite  difference  equations  and  A  has  the 
form  given  in  equation  (2-4),  With  the  cross-derivative  given  by  equation  (2-1)  E  is 
an  (N  -  1)  x  (N  -  1)  matrix  with 

E  ■  -  2a  +  2c  -a-  d/2N 

-a  +  d/2N  O 

®  NSS\^^^  -a  -  d/2N 

-a  +  d/2N  '  2a  +  2c 

and  F  is  given  in  equation  (2-5).  With  the  cross-derivative  as  given  in  equations 
(3-1) 


E  -  -  /  2a  -  2b  +  2c  -a  +  b  -  d/2N 


—  a  +  b  + 


o 
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when  b  > 0  .  The  corresponding  matrices  when  b  <  0  can  easily  be  derived. 

Equation  (3-3)  will  be  solved  iteratively  by  a  successive  line  over-relaxation 

(SLOR)  scheme,  in  which  all  values  of  $  on  a  line  y  ■  constant  are  updated  at 

the  same  time.  At  points  before  and  on  the  line  where  <J>  is  currently  being  updated 

there  is  a  choice  as  to  whether  to  use  the  values  of  d>  from  the  current  iteration,  or 

values  of  4>  found  during  the  previous  iteration.  The  first  derivative,  ♦  ,  will  be 

calculated  using  values  of  <f>  found  during  the  previous  iteration,  while  the  second 

derivatives,  4>  and  4  ,  will  be  estimated  using  values  of  $  from  the  current 

xx  yy 

iteration,  whenever  possible.  For  the  mixed  second  derivative,  4>  ,  values  of  4> 

xy 

from  the  current  iteration  will  be  used  on  the  line  previous  to  the  current  one,  while 
on  the  current  line  various  schemes,  depending  on  the  representation  of  ♦  ,  will  be 

investigated.  These  schemes  are  as  follows: 


Scheme  I 
used. 


4> 


xy 


is  given  in  equation  (2-3). 


No  values  of 


on  the  current  line  are 


Scheme  II  <t _  is  given  in  equations  (3-1).  The  values  of  4>  currently  being  calcu- 

xy 

lated  are  used  everywhere  on  the  current  line. 

Scheme  III  is  given  in  equations  (3-2) .  The  values  of  ?  current'. y  being  calcu- 

xy 

lated  are  used  everywhere  on  the  current  line. 

Scheme  IV  #  is  given  in  equations  (3-1).  The  value  of  $  currently  being  calculated 

'  *y 

is  used  at  the  point  at  which  the  derivative  is  centred.  Elsewhere  on  the  current  line 
values  of  $  from  the  previous  iteration  are  used. 

Scheme  V  $ _ is  given  in  equations  (3-2).  The  value  of  <f>  currently  being 

calculated  is  used  at  the  point  at  which  the  derivative  is  centred.  Elsewhere  on  the 

current  line  values  of  $  from  the  previous  iteration  are  used. 
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Schema  VI  $  is  given  in  equations  (3-1).  Values  of  $  from  the  previous  iteration 

~  *y  *” 

are  used  to  approximate  the  third  and  fourth  terms.  Values  of  $  currently  being 
calculated  are  used  in  the  fifth  and  sixth  terms. 


If  d  *  0  and  the  matrix,  E  ,  is  determined  from  equation  (2-4)  and  any  one 
of  equations  (2-3),  (3-1)  and  (3-2)  it  is  easily  seen  that  A  is  hermitian  and  E  is 
'$  ;  negative  definite.  Then,  if  any  of  schemes  I,  II  and  III  is  employed  it  can  be  shown 

that  the  scheme  will  be  stable,  in  the  sense  of  all  errors  eventually  decaying,  if  the 
t  relaxation  factor,  to  ,  satisfies  0<  u  <2  and  the  corresponding  matrix,  A  ,  is 

t  .8 

t  negative  definite  .  This,  however,  seems  cf  little  use,  as  this  definition  of  stability 

l  has  already  been  shown  to  be  of  limited  value  (see  section  2.4.2).  Also  d  will  often 

t  not  be  zero. 

b  :  ■ 

I  A  von  Neumann  type  analysis  can,  however,  be  used  to  give  some  indications  of 

stability  in  the  sense  of  all  the  definitions  given  in  section  1.  The  amplification 


udlLd^S^iiiilfUuH  . u-s~ i . tisb.jl.it.gi; >..k Jl.-Jn 


Case  (1)  a  and  c  of  same  ordtor  of  magnitude 
The  stability  criterion  |X|  <  1  gives 


d2  ^  4(2  -  to)2(ac  ~  b2) 
<d(4  o)l> 


Schemas  I,  11  and  111 


<  ^U(2  -  u)  -  2b]  2(ac  -  b2) 
N2  oic  [c4  -  id  -  4b] 


and  [c2  -  id  -  2b]  >  0 

Scheme  IV  (b  >  0) 


fc. 


r 

F 

i- 

$ 


■ 


d2  ^  4  [c(2  -  (p)  ±  2b]  2 (ac  -  b2) 
n2  wcTc(4  -  tjj)  +  4b]  ~ 


b  >  0  Scheme  V 


d^  _  4db(2  -  a)  <  4(2  -u)2(ac  -  b2) 
N2  N(|4  -  o>)  u) <! A  -  id) 


b  >0 


Scheme  VI  . 


The  lack  of  symmetry  in  scheme  VI  seems  to  arise  from  the  lack  of  symmetry  in  the 
representation  of  the  cross-derivative,  for  example  through  using  an  old  value  for  ,  . 

i+ 1,  j 

while  using  the  current  value  of  $.  .  .  in  equation  (3-la).  Scheme  V  is  less  restric- 

l—i  ,j 

tive  than  Schemes  I,  II  and  III,  but  Scheme  IV  is  more  restrictive  than  all  these  schemes. 

All  the  schemes  were  tested  numerically.  <p  was  initially  set  to  100.0  everywhere, 
and  the  scheme  was  said  to  have  converged  at  the  nth  iteration,  where  n  is  the  smallest 
integer  satisfying 

max  I  ~  ^  *  I  ^  1 0  7 
k 


where  i)n  is  the  vector  giving  the  estimates  of  $  after  the  nth  iteration.  The  results 

obtained  when  a  ■  c  ■  1  are  shown  in  Tabid  1  . 

For  the  case  b  -  0.9  and  d  -  0  there  seems  to  be  little  to  choose  between 
Schemes  I,  II  and  III,  but  the  other  schemes  are  noticeably  worse.  For  Scheme  IV  the 
condition  c(2  -  id)  > 2b  is  violated  and  so  the  iterative  scheme  is  divergent.  When 
b  is  0.45  it  is  predicted  that  Scheme  IV  will  be  stable  if  and  only  if  id  < 1 . 1  .  In 
fact  when  N  is  10  and  id  ■  1.1 l  the  scheme  is  convergent,  although  when  N  is 
increased  to  50,  with  id  still  1.11  the  scheme  is  divergent.  This  illustrates  the  point 
that  schemes  which  may  appear  satisfactory  for  large  step  lengths  become  less  and  less 
satisfactory  as  the  step  length  decreases. 

Attention  was  then  turned  to  d  non-zero.  The  above  analysis  suggests  that  when 
b  *  0 . 9 ,  a  ■  c  «*  1  and  id  ■  1.5 

| d/N |  <  0.225  Scheme  I 

! d/N |  <  0.439  Scheme  V 

for  stability.  Although,  as  expected,  Scheme  V  permits  larger  values  of  | d/N [  than 
Scheme  I  the  results  given  in  the  table  are  clearly  in  need  of  some  explanation! 
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To  try  to  do  so  the  case  a 
The  above  analysis  suggests 


<j  ■  id  ■  1  and  b  •  0  ,  with  N 
| d/N  |  <  1.15 


50  was  investigated. 


is  a  necessary  condition  for  stability  in  all  the  senses  of  section  I.  If,  however, 
stability  is  defined  as  occurring  if  all  errors  eventually  decay,  then  stability  will 
occur  if 

max  |  y  I  <  1 
s,t 

where  y  satisfies 

n 

2y  -  ly - j  cos  6  -  /y  cos  ■  0 

-v  4N 

with  0  *  tts/N  and  <f>  ■  irt/N  ,  s  and  t  being  positive  integers  satisfying 
1  <  s,t  <  N  -  1  .  As  N  -+■  *  this  gives 

| d/N |  <  1.82  . 


For  d/N  just  greater  than  1.15  the  iterative  scheme  works  well,  but  as  d/N  approaches 

1.82  the  initial  ri«e  in  max  |  is  so  large  that  it  becomes  unacceptable,  and 

k  x  X 

also  causes  a  large  increase  in  the  number  of  iterations  required.  Just  where  the 
iterative  scheme  becomes  unacceptable  is  not  easy  to  define,  but  it  is  clearly  here 
being  ultra-cautious  to  require  | d/N |  <1.15  ,  as  the  intermediate  results  are  not  of 
interest,  although  allowing  it  to  get  too  near  1.82  is  unacceptable. 

Case  (ii)  An  example  when  a  and  c  are  not  of  the  same  order  of  magnitude 

2  2 

It  wi  11  be  assumed  that  a0  ►  c  (**  a0  ►  b  because  b  <  ac) 


b  (d/N) 


c  «  (d0/N) 


and  thct  there  is  an  upper  limit  on  the  3ize  of  N  .  This  situation  can  arise  in  the 
soluti  m  of  the  full  potential  equations  of  fluid  motion  round  very  highly  swept 
wings  us’ng  a  non-orthogonal  grid,  in  some  regions  of  which  the  angles  between  coordinate 
lines  of  different  families  are  small. 

A  von  Neumann-type  analysis  gives  approximately, 

|d|  <  2cN  tan  | Schemes  I,  II,  III  and  VI  (3-4) 

|d|  <  2 all  tan 

for  stability,  in  all  the  senses  of  section  1.  No  stability  criterion  is  possible  for 
Scheme  IV  for  the  following  reason.  The  SLOR  method  requires  the  solution  of  equations 

i  .  .  j 


”  * 


I  + 


a  tan^  J  (2  -  w) 
N 


Scheme  V 


(3-5) 
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of  the  fora 


Tx 


(3-6) 


where  T  is  the  (N  -  1)  *  (N  -  1)  matrix  of  form 


"n-1 


with  a^  ■  c.  ■  -  a  and  ■  2(i  -  |b|  +  c). 
This  is  done  by  setting 


bl  ’ 


c. 

1 


8J  ■  bj  * 


b.  -  a. to.  . 
l  l  l-l 


ki  “  aigi-l 
b.  -  a.*.., 


(3-7) 


the  solution  being  obtained  from 

Vi  " 


%-l 


x. 

l 


8,- 


w , x,  . 

l  i+l 


(3-8) 


As  a  >  | b |  >  c  >  0  ,  there  is  nothing  to  prevent  |w.J  becoming  very  large  occasionally. 
If  |(o.|  is  very  large  t  effect  of  rounding  errors  will  not  be  negligible  and  the 
scheme  will  be  unstable.  A  sufficient  condition  to  prevent  any  SLOR  scheme  being  unstable 
in  this  way  is  to  require  the  matrix,  T  ,  to  be  diagonally  dominant,  for  this  ensures 
|uk|  <  1  for  all  i  . 

.  .  4 

The  results  obtained  when  a  ■  10  ,  b  ■  90  and  c  ■  1  are  shown  in  Table  2. 

Scheme  IV  appears  satisfactory  when  N  ■  10  and  d  ■  0  ,  but  was  found  to  be  divergent 
when  N  ■  50  and  d  ■  0  .  There  is  little  to  choose  between  Schemes  I,  II,  III  and  VI 
but  Scheme  V  is  definitely  slightly  slower  (though  not  so  markedly  as  in  the  case 
a  ■  c  *  1,  b  ■  0.9  and  d  »  0) . 

The  symmetrical  nature  of  Schemes  I,  II  and  III  and  the  unsymmetrical  nature  of 

Scheme  VI,  with  respect  to  the  sign  of  d  ,  are  illustrated,  |d|  arbitrarily  being 
4 

chosen  as  5.5  *  10  .  The  criteria  (3-4)  and  (3-5)  suggest  that,  for  stability, 

i  i  4 

when  | d |  ■  5.5  *  10  ,  in  all  the  senses  of  section  1, 


o 

00 


to  <  1.13 


Schemes  I,  II,  III  and  VI 
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2  + 


w  ^ 


1  + 


a  Can 


2  w 


4a 


T!T~  nr 

N  tan  rr 


N 


Schama  V 


00 

o 


This  implies  that  all  values  of  cu  are  permissible  with  Scheme  V  -  the  numerical  results 
support  this.  However,  the  restriction  on  the  relaxation  factor,  u  ,  does  not  seem  to 
be  quite  correct  for  the  other  schemes. 

To  try  to  explain  this,  it  should  first  be  observed  that  the  criterion  ensuring 

that  all  errors  eventually  decay  is  also  given,  approximately,  by  equation  (3-4)  for 

Schemes  I,  II,  III  and  VI.  This  suggests  that  as  the  relaxation  factor,  w  ,  approaches 

its  theoretical  upper  bound  the  constant,  of  the  second  stability  criterion, 

equation  (2-10),  becomes  too  large.  This  permits  an  unacceptably  large  initial  growth 

of  errors,  although,  for  fixed  u  ,  no  error  can  become  arbitrarily  large.  Thus  the 
i  i  4 

results  when  |d|  -  5.5  x  10  illustrate  the  possibly  unsatisfactory  nature  of  the 
second  stability  criterion.  Once  again,  just  where  the  iterative  scheme  becomes 
unsatisfactory  is  difficult  to  define. 

For  a  larger  value  of  d  -  a  value  of  1  *  10^  was  taken  -  with  Scheme  V,  the 
criterion  (3-5)  implies,  approximately 

u  <  1.21 

for  stability  in  all  the  senses  of  section  1 ,  However,  if  it  is  merely  required,  for 
stability,  that  all  errors  should  eventually  decay,  the  stability  criterion  may  be 
relaxed  to,  approximately 

8a  |a  sin2  ^  +  bJ 


w 


< 


dV 


indicating  that  all  values  of  the  relaxation  factor,  u  ,  will  be  permissible.  In  the 
numerical  work  it  was  found  that  values  of  w  somewhat  greater  than  1.21  were 
acceptable,  but  as  w  increased  the  initial  rise  in 


max 

k 


?k 


,n-l 


is  so  large  that  the  scheme  becomes  unacceptable.  Just  where  the  iterative  scheme 
becomes  unsatisfactory  is  difficult  to  define.  However,  as  the  intermediate  results  are 
not  of  interest,  it  is  clearly  too  cautious  to  require  u  <  1.21  ,  although  allowing  it 
to  get  too  large  is  unsatisfactory. 

Discussion 

The  results  of  this  section  have  illustrated  many  of  the  points  made  in  section  2, 
concerning  possible  stability  criteria.  A  von  Neumann-type  stability  analysis  has  been 
shown  to  be  of  considerable  use  -  although  it  has  shortcomings,  some  of  which  have  been 


triaeiag3!fa-iH- 


'•if- 1  t=v  -rr ; 


4 
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illustrated .  In  most  of  the  section  the  von  Neumann  criterion  haa  bean  aeen  to  be  | 

unneceeaarily  stringent,  although  in  on*  caae  it  waa  not  quit*  atringant  enough.  | 

i 

The  reaulta  also  suggest  ways  in  which  mixed  partial  derivatives  should  be  handled  J 

in  second-order  elliptic  partial  differential  equations  when  SLOR  is  being  used  to  1 

V5j 

aol/e  the  finite-difference  equations.  Concerning  the  number  of  iterations  required,  and 

regions  of  stability,  there  is  little  to  choose  between  Schemes  1,  II  and  III.  However,  j 

the  calculations  for  Scheme  I  vill  take  slightly  less  time  than  those  for  Schemes  II  and 

III.  Schemes  V  and  VI  require  more  iterations  than  Schema  I,  but  Scheme  V  is  stable  in 

circumstances  where  Schemes  I,  II  and  III  are  unstable.  The  region  of  stability  of 

Scheme  IV  is  very  much  less  than  the  region  of  stability  of  Schemes  I,  II  and  III,  so 

Scheme  IV  is  of  relatively  little  use.  The  region  of  stability  of  Scheme  VI  is  about 

the  same  as  the  regions  of  stability  of  Schemes  I,  II  and  III,  so  Scheme  VI  is  worse  q 

than  Scheme  I.  Thus  it  is  recommended  that  Scheme  I  should  be  used  whenever  possible, 


but  if  Scheme  I  is  unstable  Scheme  V  should  be  tried. 

3.2  Example  2  -  first  derivative 

In  this  example  the  manner  in  which  instabilities  can  arise  from  a  boundary 
condition  will  be  considered.  A  suitable  method  for  handling  some  of  the  first 
derivative  terms  which  arise  in  a  second-order  differential  equation  will  also  be 
indicated. 

In  the  solution  of  a  second-order  partial  dif  erential  equation  there  may  be 
regions  in  which  the  equation  reduces  to,  essentia. ly, 


with  homogeneous  boundary  conditions.  This  will  arise  through  either  the  coefficients 
of  the  other  derivatives  being  small,  or  through  derivatives  in  directions  other  than 
the  x  direction  being  small  -  as  can  happen  in  the  solution  of  the  full  potential 
equations  of  fluid  motion  round  very  highly  swept  wings  when  using  a  non-orthogonal 
grid.  In  such  a  situation,  the  manner  in  which  the  boundary  conditions  can  cause 
instabilities  is  illustrated  by  the  following  example!  solve  (3-9)  on  (0,l)  with 


0  at  x  ■  0  ,  b  ■  pt  at  x  ■  1 ,  p  >  1 


(3-10) 


Take  the  usual  central  difference  representations  of  *x  and  .  Let  Ax  -  1/N 

and  let  the  subscript  i  refer  to  the  coordinate  direction,  x  .  The  finite  difference 
equations  may  be  written  in  the  forms 


2  -1- 


V 

2  -  ■?£  I* 

Z  N  \9N 


(3-11) 


where  is  the  estimated  value  of  $(i/N)  . 
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Several  schemes  will  be  considered  to  solve  this  equation.  They  are  all  of  the 

formi 

find  vn+1  from  Hvn+I  ■  Gvn  (3-12) 

where  H  -  G  ■  A  ,  and  apply  over-'elaxation. 

Scheme  I  Here  the  second  derivative  will  be  estimated  from  the  values  of  $  currently 
being  calculated.  Hence  H  is  the  same  as  A  of  equation  (3-11)  and  G  -  0  . 

Scheme  II  Even  when  equation  (3-9)  is  a  good  approximation  to  the  full  equation,  the  full 
equation  itself  may  be  much  more  complex.  In  such  circumstances  it  is  quite  likely  that 
some  of  the  first  derivative  term  will  be  evaluated  using  values  of  $  currently  being 
calculated  while  the  rest  of  the  first  derivative  will  be  estimated  using  values  of  $ 
from  the  previous  iteration.  This  is  best  understood  by  writing  equation  (3-9)  as 


xx 


■ 

6  X 


8*, 


where  |g/2Nj  is  not  necessarily  small. 

Terms  on  the  left-hand  side  are  to  be  estimated  using  the  values  of  $  currently  being 
calculated,  while  the  term  on  the  right-hand  side  is  to  be  found  using  values  of 
from  the  previous  iteration.  H  is  now 


and  G  is 


The  method  of  matrix  inversion  given  in  equations  (3-6),  (3-7)  and  (3-8)  will  be 

n+  j 

used  to  find  y  in  equation  (3-12).  To  ensure  that  this  is  stable,  |g|  <  2N  -this 

ensures  that,  with  the  possible  exception  of  the  last  row,  H  is  diagonally  dominant. 
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Scheme  III  If  |g|  >  2N  Scheme  II  will  not  be  satisfactory.  However,  now  write 
equation  (3-9)  in  the  form 


*  +  g*  ~  N(g  —  2  N)*  ■  g*  -  N(g  -  2N)$ 

XX  X  x 

and  evaluate  terms  on  the  left-hand  side  using  the  values  of  $  currently  being  calcu¬ 
lated,  but  terms  on  the  right-hand  side  using  the  values  of  d>  found  during  the  previous 
iteration.  Then 


1  -  g,/2N 


-  1  +  g/2N 


1  +  g/2N  ^-g/N  -  1  -  g/2N 


-2  g/N  -  (1  +  g/2N)  2p/N  , 


and  G  ■  / g/N  -  2 


'g/N  -  2 


-  g/2N 


g/N  -  2  -  g/2N  2p/N 


If  g  >  2N  the  method  given  in  equations  (3-6),  (3-7)  and  (3-8)  for  finding  yn  in 
equation  (3-12)  is  stable,  because  H  is  diagonally  dominant  when  its  last  row  is 
ignored. 

The  stability  of  the  three  schemes  will  now  be  considered.  A  von-Neumann  analysis 
gives  the  following  equations  for  the  salification  factors,  X  : 


Scheme  I 


Scheme  II 


<U  -  1)  (A  +  u  -  1) 


(U  -  0“ 


(3-13) 


(U  -  D  +  -A-  (P  +  0 


x 


i 


(3-14) 


f 
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Scheme  III 

x  .  i _ <K  -  rf* _ 

(y2  +  I)  +  ^  (y2  -  2y  -  1) 


(3-15) 


where  | y  |  *  1  . 

In  none  of  these  schemes  can  |x|  be  greater  than  unity  for  any  permitted  value  of  y  , 
so  the  von-Neumann  stability  criterion  is  satisfied. 

The  boundary  condition  at  x  -  0  trivially  satisfies  the  Godunov-Ryabenkii 
criterion.  As  the  boundary  condition  at  x  ■  1  is  being  imposed  implicitly  the  finite- 
difference  equation  modelling  this  boundary  condition  is  independent  of  X  .  It  is 


2 

y 


1 


For  modes  which  decay  away  from  the  boundary  |y]  >  1  ,  hence  the  only  solution  it  is 
necessary  to  consider  is 

y 

which  is  real  and  greater  than  unity. 

it  is  found  that  X  «  1  -  ai  and  so  |x|  <  1  .  Thus  the  first  scheme  satisfies  the 
Godunov-Ryabenkii  criterion.  As  y  >  1  and  0  <  u  <  2 


Substituting  this  value  of  y  into  equation  (3-13) 


0  <  - CjLjLJh -  <  2 

(y  ~  i) +  2N  (y + 1) 

and  so,  from  equation  (3-14),  Scheme  II  also  satisfies  the  Godunov-Ryabenkii  criterion. 
However,  in  equation  (3-15),  |x|  will  exceed  unity  if 

2  (y2  -  2y  -  1)  +  2(y2  +  1)  <  u(y  -  l)2 

which,  for  large  g/2N  ,  gives,  approximately  y  <  1  +  /I  which  implies  that  p/N  <  1  . 

The  problem  given  in  equations  (3-9)  and  (3-10)  was  solved  numerically  in  the 
case  N  =*  50  ,  with  a  relaxation  factor  of  1.6.  In  Scheme  II  g  was  taken  to  be  100  and 
in  Scheme  III  5000,  The  results  are  shown  in  Table  3,  convergence  being  defined  as  in 
section  3.1.  They  are  in  good  agreement  with  the  above  theory. 
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Discussion 


The  results  of  this  section  have  illustrated  some  of  the  points  made  in  section  2, 
about  possible  stability  criteria.  As  the  von  Neumann  criterion  totally  ignores  the 
boundary  conditions  it  is  possible  that  even  if  the  von  Neumann  criterion  is  satisfied 
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the  scheme  may  still  be  unstable  as  a  result  of  the  particular  boundary  conditions 
used. 

The  results  also  suggest  ways  in  which  first  derivatives  should  be  handled  in 
second  order  partial  differential  equations.  It  seems  safest  always  to  use  values  of 
the  dependent  variable  from  the  previous  iteration.  If  this  is  not  done,  too  much  of 
the  first  derivative  may  be  evaluated  using  the  current  values  of  $  (Scheme  II  is 
unsatisfactory  if  |g|  >  2N) .  If  this  is  avoided  (as  in  Scheme  III)  the  boundary 
conditions  may  introduce  instabilities. 

4  CONCLUSIONS 

The  question  of  the  stability  of  iterative  schemes  has  been  discussed  with  the 
aid  of  numerical  examples.  It  has  been  shown  that  there  is  in  general  no  satisfactory 
theoretical  definition  of  stability.  This  necessarily  means  that  in  general  there  can 
be  no  satisfactory  criterion  for  stability. 

In  many  schemes  linear  equations  of  the  form  Tx  °  k  must  be  solved.  The  first 
requirement  for  stability  is  that  the  method  employed  to  solve  these  equations  must  be 
stable.  The  usual  method,  if  T  is  tridifloOnal,  of  solving  these  equations,  is  given  in 
equations  (3-6),  (3-7)  and  (3-8).  This  method  is  stable  if  T  is  diagonally  dominant. 
However,  the  results  of  section  3.2  show  that  T  being  diagonally  dominant  is  not 
always  necessary. 

If  a  stable  method  of  solving  equations  of  the  form  Tx  =  k  is  used,  a  possible 
fairly  simple  criterion  for  stability,  is  that  of  von  Neumann.  The  examples  illustrate 
that  this  criterion  can  be  unnecessarily  severe  for  some  problems,  and  not  sufficiently 
severe  for  others. 

The  first  example,  concerned  with  the  representation  of  a  mixed  second  derivative 
(section  3.1),  shows  that  if,  in  practice,  there  is  a  lower  bound  on  the  size  of  the 
step  length  taken  then  the  criterion  may  be  too  severe.  However,  in  this  example,  it 
was  only  the  precise  borderlines  between  practical  stability  and  instability  that  were 
incorrectly  predicted,  all  relative  trends  within  and  between  schemes  being  correctly 
predicted. 

The  second  example,  concerned  with  a  first  derivative  (section  3.2),  shows  that 
instabilities  may  arise  through  the  boundary  conditions,  so  that  the  criterion  (which 
ignores  the  boundary  conditions)  may  not  be  severe  enough.  Thus  the  criterion  indicates 
that  stability  may  be  obtained  if  the  boundary  conditions  are  suitable.  As  it  ignores 
the  boundary  conditions  it  cannot  tell  which  boundary  conditions  will  give  stability,  and 
which  not. 

Thus,  while  the  von-Neumann  criterion  may  be  of  considerable  help  in  choosing  a 
suitable  numerical  scheme,  it  must  be  applied  with  great  care  and  its  shortcomings  kept 
in  mind . 

The  examples  of  section  3  suggest  how  mixed  and  first  derivatives  should  be  handled 
in  SLOR  (successive  line  over— relaxation)  schemes.  If  the  mixed  derivative  is  required 
at  the  point  x  ■  iAx  ,  y  ■  jAy  then  the  values  of  the  variable  at  the  points 
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x  ■  (i  ±  l)Ax  ,  y  -  (j  ±  l)Ay  should  be  used  if  possible.  However,  if  stability 
cannot  be  obtained  with  this  scheme,  a  more  complex  one  (Scheme  V  of  section  3.1)  may 
give  stability.  When  evaluation  of  a  first  derivative  is  required,  care  should  be  taken 
if  some  first  derivative  terms,  on  the  line  currently  being  updated,  are  evaluated  using 
values  of  <j>  currently  being  calculated,  while  others  are  evaluated  using  values  oi:  <f> 
from  the  previous  iteration.  Stability  is  more  likely  to  be  obtained  if  first  derivatives 
are  always  estimated  using  values  of  4>  from  the  previous  iterations. 
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Table  1 

Nuaber  of  Iterations  required  to  solve  ♦  +  2b*  +  *  +  d*  ■  0 

_  xx  xy  yy  x 


N 

b 

d 

Scheme 

Relaxation  factor, 

0) 

Number  «'f 
iterations 

10 

0.9 

0 

I 

1.45  (optimum) 

31 

ti 

tl 

tt 

II 

1.46  (optimum) 

7 

ti 

tl 

tl 

III 

1.45  (optimum) 

30 

»t 

tl 

tt 

IV 

1.46 

D 

it 

II 

tt 

V 

1.54  (optimum) 

158 

tt 

tl 

It 

tt 

1.45 

175 

it 

tl 

It 

VI 

1.45  <  to  <  1.60 

70-73 

t» 

0.45 

II 

IV 

1.10 

157 

tt 

tt 

tl 

tt 

1.11 

182 

•t 

tt 

tt 

11 

1.20 

D 

50 

tt 

tt 

II 

1.11 

tt 

10 

0.9 

6.0 

I 

1.5 

218 

tl 

tl 

7.5 

It 

tt 

D 

tt 

tt 

tt 

V 

II 

91 

50 

0 

55.0 

1.0 

130 

It 

It 

60.0 

o 

1 

II 

145 

II 

70.0 

43 

It 

203 

ti 

" 

80.0 

eft 

<d 

U 

tl 

363 

tt 

tt 

85.0 

8 

5 

tl 

617 

tt 

tl 

90.0 

fH 

0) 

t 

•H 

II 

2355 

It 

It 

95.0 

II 

D 

Key:  D  denotes  divergence 
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Table  2 


4 

Number  of  iteration*  required  to  solve  10  ♦ 


Scheme 


Relaxation  factor. 


1.C0  (optimum) 


1.08  (optimum) 


1  <  1.95 


1  <1  <1.95 


Number  of 
iterations 


D  denotes  divergence 

C  denotes  convergence,  but  case  not  run  to  full  convergence 


Table  3 


Number  of  iterations  required  to  solve  4^ 


with  4-0  at  x 


0  and  4 
_  x 


p4  at  x 


-  J 


■ 

Scheme 

Number  of 
iterations 

25/? 

I 

43 

II 

II 

C 

II 

III 

D 

45 

II 

C 

II 

III 

D 

51 

II 

C 

II 

III 

C 

Key:  D  denotes  divergence 

C  denotes  convergence,  but  case 
not  run  to  full  convergence 
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